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We constructed two sets of twist-3 correlation functions that are responsible for generating the 
novel single transverse-spin asymmetry in the QCD collinear factorization approach. We derive 
evolution equations for these universal three-parton correlation functions. We calculate evolution 
kernels relevant to the gluonic pole contributions to the asymmetry at the order of as- We find 
that all evolution kernels are infrared safe as they should be and have a lot in common to the 
DGLAP evolution kernels of unpolarized parton distributions. By solving the evolution equations, 
we explicitly demonstrate the factorization scale dependence of these twist-3 correlation functions. 
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I. INTRODUCTION 

Large single transverse-spin asymmetries (SSAs), = {a{sT) — cr(— .st))/(o'(st) +o'(^'St)); defined as the ratio of 

I the difference and the sum of the cross sections when the spin st is flipped^ave been consistently observed in various 
experiments involving one polarized hadron at different collision energies [U, Q . From the parity and time-reversal 
invariance of the strong interaction dynamics, the measured large asymmetries in high energy collisions should be 
directly connected to the transverse motion of partons inside a polarized hadron. Experimental data on SSAs provide 
excellent opportunities to probe QCD dynamics beyond what have been explored by the very successful leading power 
QCD collinear factorization formalism [j, 0, @|. 

For high energy scattering cross sections with a large momentum transfer, Q ^ Aqcd, QCD collinear factorization 
approach is more relevant, and SSAs of these cross sections could be generated by twist-3 multi-parton correlation 
functions (tI. [sI. [qI. [lOl. . With the parton's transverse momentum integrated, these multi-parton correlation functions 

' represent a net (or integrated) spin dependence of parton's transverse motion inside a polarized hadron. On the 
other hand, SSAs of cross sections with two different momentum transfer scales, Qi ^ Q2 Aqcd, could be 
expressed in terms of the transverse momentum dependent (TMD) parton distributions, which directly probe the 
spin dependence of parton's transverse motion at the momentum scale Q2 while the larger scale Qi defines the hard 
collision [l^. ITsl. IT^. iTsl. flGj . These two approaches each have their own kinematic domain of validity and were shown 

] to be consistent with each other in the kinematic regime where they both apply ^17'|. Both approaches have been 
applied extensively in phenomenological studies [13, [H [13, [lH, [13, [23,, .24, .25, 2(3, 21, ,2^, 29, M, 31, M , and have 

I had initial successes [3l,[33]. However, all existing perturbative calculations are at the leading order (LO) in strong 
coupling constant, as{n), and have a strong dependence on the choice of the renormalization scale /x as well as the 
factorization scale fip, while the physically observed SSAs should not depend on the choice of the renormalization 
and/or the factorization scale. The strong dependence on the choice of renormalization and factorization scale is an 

' artifact of the lowest order perturbative calculation, and a significant cancellation of the scale dependence between the 
leading and the next-to-leading order (NLO) contribution is expected from the QCD factorization theorem [1, [1, [1, [1] . 
In order to test QCD dynamics for SSAs, it is necessary to calculate the evolution (or the scale dependence) of the 
universal long-distance distributions and to evaluate the perturbative short-distance contribution beyond the lowest 
order. 

Since the 80's, a tremendous effort has been devoted to derive evolution equations of twist-3 correlation functions 
that contribute to the structure function 52, which can be extracted from measurements of the double-spin asym- 
metry for cross sections of inclusive deep inelastic scattering (DIS) between a longitudinally polarized lepton and 
a transversely polarized hadron [51, [sl] . Evolution equations for the chiral odd twist-3 correlation functions of an 
unpolarized hadron and e, which are relevant to the SSAs in connection with the twist-2 transversity distribution, 
were also derived by several authors [STf . Evolution of the moments of TMD parton distributions was discussed in 
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Ref. [38|. However, a systematic study of the factorization scale dependence for the set of twist-3 correlation functions 
that are responsible for the gluonic and fermionic pole contributions to the SSAs 0, H, EH is not available in the 
literature. It is important to note that this set of twist-3 correlation functions has a close connection to the TMD 
parton distributions and does not contribute to the inclusive structure function 92 • In most existing phenomeno- 
logical studies of SSAs, the scale dependence of these correlation functions was often assumed to be the same as that 
of the unpolarized parton distribution functions (PDFs) [1, [HI, HH, US, ■ 

In this paper, we construct two sets of twist-3 correlation functions that are responsible for generating the SSAs in 
the QCD collinear factorization approach, and derive a closed set of evolution equations for these universal twist-3 
multi-parton correlation functions. We then focus on the evolution of two types of twist-3 multi-parton correlation 
functions that give the leading gluonic contribution to SSAs 0, The first type covers the quark-gluon correlation 
functions defined as 0, [1| , 

^,{y^)\P,ST). (1) 

There is one quark-gluon correlation function, Tg p, for each quark (anti-quark) flavor q (q). The second type covers 
the correlation functions of three active gluons [23, [28| , 

F"+(yr)|P,.T). (2) 

There are two independent tri-gluon correlation functions, Tqp{x,x) and Tq'p(x,x), because of the fact that the 
color of the three gluon field strengths in Eq. ^ can be neutralized by contracting with either the antisymmetric ifabc 
or the symmetric dabc tensors with color indices, a, 6, and c p?!. [28j. In above equations, the proper gauge links that 
ensure the gauge invariance of these correlation functions have been suppressed [8| . The jip \s the factorization scale 
and e^Tcmn _ e^'^'^'^ STpni^^h^ with the light-cone vectors, ~ {n^ ^ut) ~ (0,1, Ot) and — (1,0, Ot), which 
project out the light-cone components of any four- vector asV ■ n — and V ■ n — V" . In Eqs. ([T]) and (l2|), the 
subscript "i^" refers to the fact that it is a field strength operator F^'^ (rather than a covariant derivative operator 
D„) in the square brackets that represents the middle active parton. In this paper, we also calculate the evolution 
kernels for these two types of correlation functions at the first non-trivial order in as(/i) and study the factorization 
scale /x^-dependence of these correlation functions. 

There could be many different approaches to derive the evolution equations for the factorization scale dependence 
of these twist-3 correlation functions. Since these correlation functions are universal, the evolution equations should 
not depend on how they were derived. In this paper, we derive the evolution equations from the Feynman diagram 
representation of these correlation functions [39|. We first introduce the Feynman diagram representation for these 
twist-3 correlation functions that are relevant to the SSAs. From the operator definition of the twist-3 correlation 
functions, we then derive the cut vertices in momentum space to explicitly connect these correlation functions to 
Feynman diagrams [40| . Following the technique introduced in Ref. [3^ . we derive the evolution equations in two 
steps. First, we factorize, in terms of QCD collinear factorization approach [l,|4l|,|43, the perturbative modification to 
the twist-3 correlation functions into a convolution of the short-distance evolution kernels with the twist-3 correlation 
functions. Then, we calculate corresponding evolution kernels in the light-cone gauge. We also provide the prescription 
to calculate the evolution kernels in a covariant gauge which should give the same results. 

The existence of evolution equations of these correlation functions is an immediate consequence of the QCD collinear 
factorization of the single transverse-spin dependent cross sections [H, 13] . The general form for a factorized hadronic 
inclusive cross section with one hadron transversely polarized and a large momentum transfer Q may be written as 
ii 

a(Q, st) = Ho<E>f2®f2 + (1/Q) H^®h®h + 0(1/Q') • (3) 

The iJo and Hi are perturbatively calculable coefficient functions expanded in power series of a^, and /„ are nonper- 
turbative matrix elements of the products of fields on the light cone and are often loosely referred as "twist-n" parton 
distribution or correlation functions. In Eq. ([31), the "C>5" represents the convolution over partons' momentum frac- 
tions. The first term in the right-hand-side (RHS) of Eq. which is often referred as the leading twist term, does not 
contribute to the single transverse-spin dependent cross section, defined as A(t((5, st) = [c((5, st) — o'iQ, — •sr)]/2, be- 
cause of the parity and time-reversal invariance of the strong interaction. Therefore, single transverse-spin dependent 
cross section directly probe the twist-3 parton distribution or correlation functions as [3, 



Tq^Fix,X,fJ.F) = / ^TT" 
ZTT 



i;P+ 



Vl 



xP+ 



{P,st\F\{0) 



dy. 



— STO'nn 



AaiQ, st) = il/Q)HiiQ/^iF, a,) ^(Mi.) ® ^(AiF) + 0(1/Q') , 



(4) 
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where the summation over parton flavors and the dependence on renormaUzation scale have been suppressed. Since 
the physically measured cross section is independent of the choice of renormalization and factorization scale, the 
factorization scale dependence in the RHS of the factorized formula in Eq. (jlj should be cancelled between the iip- 
dependence of the short-distance coefficient functions and the /ii?-dependence of the parton distribution and correlation 
functions. The /.t^ dependence of the normal twist-2 PDFs satisfies the DGLAP evolution equation [ist . 

= (5) 

where the parton flavor dependence has been suppressed and P2 is the twist-2 evolution kernel, which can be calculated 
perturbatively and expressed in a power series of a^. From dAa{Q, ST)/dln{jj,p) ~ 0, we derive the leading order 
generic evolution equation for /a as 

' :h=(7n3-^H[^'-P^'')^h^ (6) 



by applying the factorized formula in Eq. Q on parton states and expanding it to the first non-trivial power of as- 
In deriving Eq. ([6]), we divided out the leading order coefficient function, h[^\ Equation ([6]) clearly shows that the 
evolution equation is a consequence of QCD factorization and indicates that every perturbatively factorizable single 
transverse-spin dependent cross section could be used to derive the evolution kernels of twist-3 correlation functions. 
For example, the order of evolution kernels could be obtained by calculating h[^\ the one-loop corrections to 
the short-distance partonic hard part of single transverse-spin dependent Drell-Yan cross section [4y|. The evolu- 
tion kernels derived in this way should be the same as what we derived here directly from the Feynman diagram 
representation. 

The quark-giuon and tri-gluon correlation functions in Eqs. ([T]) and ([2|) give the leadin g so ft gluonic pole contribution 
to the transverse-spin dependent cross section with a single hard scale, Acr{Q, st) [sI. [2811331. However, transverse-spin 
dependent cross sections with more than one physically observed hard scales could get additional hard pole contribution 
which is proportional to the off-diagonal part of the twist-3 correlation functions, Tq p{x, x' , ^p) and Tq p\x, x' , ^p) 
where x' is not necessarily equal to x [TtI . [2^ . In addition to the gluonic pole, the SSAs or the transverse-spin 
dependent cross sections could obtain contributions from the fermionic pole of the partonic hard scattering 0, [§1 ■ The 
leading fermionic pole contribution is generated by not only the off-diagonal part of the correlation functions Tq^p 

and Tqp^ but also a new set of twist-3 correlation functions that have a vanishing diagonal contribution [1, [l^, [IH ■ 
In order to describe the phenomenon of SSAs for observables with more than one hard scale and evaluate the full 
perturbative contribution to SSAs beyond the lowest order in as, it is necessary to study both the diagonal and 
off-diagonal twist-3 correlation functions that can generate the SSAs. 

The rest of this paper is organized as follows. In the next section, we construct two sets of twist-3 correlation 
functions that can generate the SSAs. In Sec. lIIIi we introduce the Feynman diagram representation for these twist-3 
correlation functions. To connect the Feynman diagrams to the specific twist-3 correlation functions, we derive the cut 
vertices from the operator definition of these twist-3 correlation functions in momentum space. From the perturbative 
modification to the correlation functions, we derive the evolution equations. In Sec. IIVI we calculate all evolution 
kernels at the order of for the evolution equations of the quark-gluon and tri-gluon correlation functions defined 
in Eqs. ([T]) and ([2]). In Sec.|Vl we discuss the scale dependence of these correlation functions by solving the evolution 
equations. Finally, we give our conclusions and a brief discussion of the impact of the calculated scale dependence of 
the correlation functions in Sec. IVII 



II. TWIST-3 CORRELATION FUNCTIONS RELEVANT TO SSAS 



In this section, we construct two sets of twist-3 correlation functions that are responsible forgenerating the gluonic 
and fermionic pole contributions to the SSAs in the QCD coUinear factorization approach 0, @|- 

We first introduce two twist-3 correlation functions by generalizing the definition of the diagonal functions in Eqs. ([1]) 
and ©, 

f,,pAx,x + X2,fiF,ST) ^ I ^^^l^e^'''^y^e"^P^y^{P,ST\%iO)^[F,+ {y^)] i:,{y^)\P, sr) , (7) 

and 
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where the subscript "i^" again indicates that a field strength operator (not a covariant derivative operator [1, is 
inserted in the middle of the bi-local operator that defines the twist-2 spin-averaged quark (q) or giuon (G) distribution 
function. The reality property of these two functions can be expressed as , 



T'GF\{^-^X + X2,tJ.F,ST)* = f^^p\{x + X2,X,^lF,ST)- 



(9) 



That is, the real part of these two functions arc symmetric in the exchange of x and x -\- X2, while the imaginary part 
is antisymmetric. Similarly, from the parity and time-reversal invariancc, wc find [lol | 



Tq,F,o{x,X + X2,^lF.ST) = -Tq^F,aix + X2,X,^1f,~St) : 



G.Fcy\^ + ^^i■•^^^^F■,-■ST) 



(10) 



That is, these two functions are antisymmetric when the transverse spin vector st reverses its direction. 

From the definition in Eq. ([7]) and the symmetry properties in Eqs. (j9l) and (jlOp . we construct a twist-3 quark-gluon 
correlation function that is relevant to the SSA as follows, 



r,^F{x.x + X2,tiF) = e^--"« l ^ 



Tq^F,a{x, X + X2, /iF, St) ~ Tq^F,a{x, X + X2, /iF, -St) 
%j^F,<t{x, X + X2,pLF, St) + Tq^F,a{x + ^2, X, ^1f,St 
Tq^F{x,X + X2,^J.F, St) + Tq.pix + X2,X,flF, St) 
Tq^F{x,X + X2,HF,St) 



1 
2 

Re 



(11) 



where the spin-dependent twist-3 quark-gluon correlation function is defined as 



Tq^F{.X,X + X2.,tJ,F,ST) = 



Tq^F{x + X2,x, ^iF-, -St) ■ 



K^{y2)] MyT)\P,ST), 



(12) 



As shown in Eq. pip , the twist-3 quark-gluon correlation function Tq,F{x^ x + X2, (J-f) is real and symmetric when the 
active momentum fraction x exchanges with x + X2. 

Similarly, we can construct the tri-gluon correlation function relevant to the SSA as, 



^sr.nn 1 \t^M{x, x + X2,^iF,ST)~ r^^f^ix, X + X2,^iF. ' St) 
'^G.Fm^^' X + X2,fJ.F,ST) + Tq p '^{x + X2,X, fJ,F, St) 

if 4), 



1 

2 

Re 



Tqp\x,X + X2,^J.F, st) + Tqp'{x + X2,X, ^F, St) 
T(^p'(X,X + X2,fJ-F,ST) 



(13) 



where the spin-dependent twist-3 tri-gluon correlation function is defined as 



T^■^p'{x,X + X2,^iF,ST) 



(27r)2 



^Gf \^ + ^2, X, ^iF, -St) ■ 



— {P,st\F+''{0) [e^-'^"«i^/(2/2-)] F+\yi)\P,ST){-g,>.) 

(14) 



The tri-gluon correlation function T^'^p^ {x,x + X2, fJ-p) is also real and symmetric in the exchange of x and x + X2- 

The newly defined twist-3 correlation functions, Tq,F{x,x + X2,fJ.p) and Tqp\x,x -f X2,fJ.p), are related to the 
diagonal correlation functions in Eqs. (P) and ^ as 



Tq^F{x,X, ^J,f) = J dX2 [2TrS{x2)] Tq^Fix.X + X2, ^J.F) , 

TaJpi^'^'l^F) = J dx2 [2ttS{x2)] (^^^ T^-'^p\x,x + X2, fip) ■ 



(15) 
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Notice that the new tri-gluon correlation function, Ta^p{x,x' , fip), is symmetric in the exchange of x and x', while a 
direct generalization of the diagonal tri-gluon correlation function in Eq. ([2]), Tg^f{x,x' , ^p) = Tg p(x,x',^p)/x is 
not symmetric in exchanging x and x' . 

In addition to the giuonic pole contribution, the SSA could also be generated by the fermionic pole of partonic hard 
scattering 0, The fermionic pole contribution at twist-3 is proportional to the off-diagonal part of the correlation 
fimctions Tq^p and Tcp, as well as a new set of twist-3 correlation functions which vanishes when X2 = [sl. [lol. ITll| . 
To construct this new set of twist-3 correlation functions, we introduce two new twist-3 correlation functions. 



{x,X + X2,tlF,ST) 



and 



r' 



AG,F,a 



{x, X + X2,HF, St) 



_ f dvi dy2 



P+ 



{P,st\F+p{0) [tF,+ {y^)] F+\y^)\P,ST) (*e±pA) ,(17) 



where the antisymmetric tensor e±px = e'j^ = — e''^"" and subscript "Ag" and "AG" indicate that the field strength 
operator is inserted in the middle of the bi-local field operators that define the twist-2 quark helicity distribution Aq 
and the gluon helicity distribution AG, respectively. Similar to Eq. ([9]) , the reality property of these two new twist-3 
correlation functions can be expressed as. 



TAq,P,cr{x,X + X2,fJ.p,ST)* = -T^^q,F,<j{x + X2, X , ^ip , St) , 
'^Acfpai^^^ + ^-i^f^F.ST)* = -fl^Gl^{x + X2,X,Hp,ST)- 



(18) 



That is, the real part of these two new functions are antisymmetric in the exchange of x and x + X2, while the 
imaginary part is symmetric. This reality property is different from that of the functions X;,f,(t and Tcp.cr- Similarly, 
from the parity and time-reversal invariance, we find. 



TAq.F.ai^^^ + ^-i^t^F^ST) T/:^q,F,a{x + X2, X, ^Ip , - St) ■, 

'^AG%,ai^'^ + ^2^^^F^^T) = fA^%„{x + X2,X,fJ.F,-ST)- 



(19) 



That is, these two functions are symmetric when the transverse spin vector st reverses its direction. 

From the definition of these new correlation functions in Eqs. p6p and (jl7p and their properties in Eqs. psp and 
(|19p . we construct the second set of twist-3 quark- gluon and tri-gluon correlation functions that could also contribute 
to the SSAs. The new quark-gluon correlation function is defined as, 

St T'Aq,p,a{x, X + X2, /iF, Sp) - TAq.F.ai^, X + X2,HF, -Sp) 



TAq,F{x,X + X2, ^J-p) 



1 r 



TAq,F,a{x, X + X2,/iF, Sp) - TAq,F,a{x + X2,X,flp, Sp) 
T'Aq,F{x,X + X2,HP,Sp) - TAq.pix + X2,X,fiF,ST) 
TAq,p{x,X + X2,HF,Sp) 

where the spin-dependent new twist-3 quark-gluon correlation function is defined as 



1 

2 

Re 



(20) 



TAq.,F{x, X + X2, ^J.p,Sp) 



djX^ ^xP^y^ ^x^P^y- 7-7 



[is^pF+{y^)]^,q{y^)\P,sp) 



-TAq,p{x + X2,X, flp, -Sp) 



(21) 



which was also discussed in Ref. [llj. As shown in Eq. ((20)) . this new twist-3 quark-gluon correlation function 
TAq,p{x, X + X2, fJ.p) that is relevant to the SSA is also real, but, is antisymmetric in the exchange of x and x + X2. 
Similarly, the new tri-gluon correlation function is defined as. 



n-if.d) 
^AG,F 



{x,X + X2,Mf) 



'^AG%,ai^^ ^ + ^2, MF, Sp) - Tag'},A^^ ^ + a;2, MF, - Sp) 
'^AG%,ai^' ^ + ^2, MF, Sp) - f^Gjp^^ix + X2,X, [Ip , Sp) 



if A 



1 

2 

Re 



'^AG%i^^ ^ + 3^2, MF, Sp) - T^q'p{x + X2, X, HF,Sp) 



'^AG%i^^ X + X2,tJ-F,Sp) 



(22) 
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where the spin-dependent new twist-3 tri-gluon correlation function is defined as 



(2^)2 



-^T^G%{x + X2.,X,^iF,-ST) (23) 



From Eq. (j22p . it is clear that the new twist-3 tri-gluon correlation function T^q''Ip(x,x + X2,^j.f) is also real, but, 
antisymmetric in the exchange of x and x + X2. Consequently, the diagonal part of these two new correlation functions 
vanishes, 



TAq.Fix,X, Hf) = J dX2 [2t: S{X2)] TAq,F{x,X + X2, fJ-F) = , 



T^J^^^p{x,x,fiF) = y dx2 [2^(5(X2)] Ti^'''(x,a; + X2,/iF) = 0. (24) 

That is, this set of twist-3 correlation functions does not directly generate soft gluonic pole contribution to the SSAs 

To complete this subsection, we summarize the key properties of these twist-3 correlation functions that are re- 
sponsible for generating the SSAs from the unpinched gluonic and fermionic poles of partonic scattering in the QCD 
coUinear factorization approach. From their operator structure, these correlation functions can be grouped into two 
sets. One set is for the Tq^F and Tqp\ and the other includes Taq.f and T^q\,. The operators for the first set of 

correlation functions, Tq^F and Tqp \ are constructed from the bi-local operators that define the twist-2 spin-averaged 
PDFs with an insertion of the following operator, 



27r ^ " ^^'^ / 27r 



e— - le-"-^- (y^-)] ^^ / ^e"^^ st, ^^.+ (2/2^)] ; (25) 



and the operators for the second set of correlation functions, Taq^f and T^^lp, are constructed from the bi-local 
operators that define the twist-2 spin- dependent parton helicity distributions with an insertion of a slightly different 
operator. 



F„^{y2)\ ■ (26) 



27r 

The it'^^ in Eq. (j25p takes care of the parity invariance of the spin asymmetry for the first set of correlation functions, 
while the same property was taken care of naturally by the 7^ or ie±px in the operator definition of the spin-dependent 
helicity distributions. The extra "i" in both Eq. and Eq. (pS)) provides the necessary phase for the SSAs and 
is a result of taking the contribution from the gluonic or fermionic pole of partonic scattering Q. It is this phase 
that both sets of twist-3 correlation functions relevant to the SSAs do not contribute to the long-distance correlation 
functions extracted from any measurement of parity conserving double-spin asymmetries. For example, none of these 
twist-3 correlation functions, 7^,f, 'Fqp\ TAq,F, and 'F^q'I^, directly contributes to the DIS structure function g2. 



III. EVOLUTION EQUATIONS 

In this section we introduce the Feynman diagram representation of the twist-3 quark-gluon and tri-gluon correlation 
functions defined in the last section. We derive cut vertices in momentum space to connect the Feynman diagrams 
to specific twist-3 correlation functions [40j . From the Feynman diagram representation, we derive the evolution 
equations for the scale dependence of these twist-3 correlation functions. 



A. Feynman diagram representation and cut vertex 

In QCD coUinear factorization approach to SSAs, the twist-3 three-parton correlation functions measure the net 
effect of the quantum interference between two scattering amplitudes of the transversely polarized hadron: one with 
single active parton and the other with two active partons, participating in the short-distance hard scattering Like 
the normal PDFs, the quark-gluon and tri-gluon correlation functions, as defined in last section, could be represented 
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FIG. 1: Feynman diagrams that contribute to the twist-3 quark-gluon (a) and tri-gluon (b) correlation functions. a,l3,^ and 
a, b, c are Lorentz and color indices of gluon field operators, respectively. 



by the cut forward scattering diagrams as sketched in Figs.[TJa) and[ljb), respectively. The cut represents a particular 
final-state. The Feynman diagrams in Fig. [1] should include all possible cuts to sum over all possible final states. We 
suppress the explicit cuts for the diagrams in Fig. [T] since the matrix element of the three- parton correlation functions 
with the middle gluon field strength in the left side of the cut is equal to the matrix element with the gluon field 
strength in the right side of the final-state cut. This is because the field operators of hadronic matrix elements commute 
on the light-cone [1,1131 . Because of the odd number of active fields defining the twist-3 correlation functions, unlike 
the normal PDFs, these correlation functions do not have the probability interpretation. 

As discussed in the last section, one set of twist-3 correlation functions is expressed in terms of a sum of two spin- 
dependent twist-3 correlation functions, as in Eqs. (|lip and (|13[) . and the other by a difference of two spin-dependent 
twist-3 correlation functions, as in Eqs. (|20p and (j22p . These spin-dependent twist-3 correlation functions are given 
in terms of explicit matrix elements of the transverse-spin dependent hadronic state and could be represented by 
Feynman diagrams. However, since all gluon lines in Feynman diagrams arc connected to gluon fields, calculating 
the Feynman diagrams in Fig. [1] does not immediately give the twist-3 correlation functions whose gluonic degree of 
freedom is represented by the field strength, not the gluon field, A'-'-. Therefore, in order to fully define the 

Feynman diagram representation of the spin-dependent twist-3 correlation functions, we need to derive the cut vertex 
[40t to connect the operator definition of the spin-dependent twist-3 correlation functions to the cut forward scattering 
Feynman diagrams in Fig. [TJ With different cut vertices, the same diagrams in Fig. [l]can represent both sets of the 
spin-dependent twist-3 correlation functions. 

The cut vertex can be derived by rewriting the operator definition of the correlation functions in terms of hadronic 
matrix elements of quark and gluon operators in momentum space [40| . 

To derive the cut vertex to connect the spin-dependent quark-gluon correlation function in Eq. (|12p to the Feynman 
diagram in Fig. [TJa) , we reexpress the operator definition of the correlation function in Eq. (|12p in terms of hadronic 
matrix elements of quark and gluon operators in momentum space and find. 



Tg^Fix,X + X2,fJ-F,ST) 



;27r)4 (27r)4 
7+ „ / k+ 



2P+ 



{P,ST\^g,,i-k-k2) 



P+ 



X2 5 X2 



kt 
P+ J 



k2an^ 



xA^'^(fc2)^,,j(fc)|P,ST), 



(27) 



where the fermionic color contraction factor Cq is given by 



(28) 



with quark and gluon color indices, i, j = 1,2,3 = Nc and c = 1,2, ...,8 = — respectively, and are the generators 
of the fundamental representation of SU(3) color. The field operators listed with "~" in Eq. ^7} represent the 

momentum space field operators of those in Eq. (fT2|) . The matrix element, {P, st\ ipq i{~k — k2)A^-''^{k2) i^qjik) \P, st), 
in Eq. (|27p can be represented by the Feynman diagram in Fig. [Ija). By comparing the definition of quark-gluon 
correlation function in Eq. ()27p and the Feynman diagram in Fig. [Ija), it is clear that we can derive the quark-gluon 
correlation function, 7^,f, from the Feynman diagram in Fig. [ija) by contracting the quark and gluon lines with the 
expression in the square brackets and the color contraction factor (t'^)ij, plus the integration over the loop momenta 
in Eq. ((27p . The expression in the square brackets plus the color contraction factor Cq defines the cut vertex that 
connects the Feynman diagram in Fig.[lja) to the quark-gluon correlation function Tq p in Eq. (|12p . 



V, 



-1+ ( k+ 

2P^^V-P^ 



') X2 S 



X2 



^2 
P+ 



k2crn^ 



c„ 



(29) 
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Similarly, we can rewrite the tri-gluon correlation function in Eq. ([M]) as 



(27r)4 (27r)' 



(P, st\ Af'-^-k - k2)A^-'{k2) A"^"(fc) \P, st) 



x{x + X2) -gap + 



(fc + k2)anp kfjUg ^ k ■ {k + k2) Ugnp 
{k + k2)+ ~k+ k+ (fc + k2)+ 



6 X 



fc+ 
P+ 



X X2S { X2 - ^ 



k2a rLf^ 
^2 



U,d) 
albca ' 



where the gluonic color contraction factor Cg is given by 

(C,)i2 = z/fcca = (^'=)6a , and {Cg)fl^d,,a.. 



(30) 



(31) 



where are the generators of adjoint representation of SU(3) color. In Eq. ([50]) . the matrix element (P, St \ Al^'^[—k — 
^2) A^''^{k2) A"'"(fc) |P, St) can be represented by the Feynman diagram in Fig.[IJb). Similar to the situation of quark- 
gluon correlation, the expression in the square brackets plus the color contraction factor C^^^'^'' defines the cut vertex 



for calculating the tri-gluon correlation function T^Q p' from the diagram in Fig.[IJb), 



V, 



G,F 



x{x + X2) -gap + 



{k + k2)anp kfjUa k ■ {k -\- k2) Tlanp 

{k + k2)+ ^ ~k+ fc+ (fc + fc2)+ 



S X- 



P^ 



xx26{x2-j^] (ze^-^' 



(Ad) 



(32) 



Similarly, by rewriting the operator definitions of the second set of spin-dependent twist-3 correlation functions in 
terms of quark and gluon field operators in momentum space, we derive the following cut vertices. 



7 7 



fc 



VAg,F = H^-—] (-St) X2S[x2-^ 



P 



P+ 



k2a1^fj 

K2 



(33) 



for connecting the same Feynman diagram in Fig. [T](a) to the second set quark-gluon correlation function 7a<j,f in 
Eq. ((2T|) . and 



Vag,f = x{x + X2) iie±p\) 
XX2S [ X2- 



(fc + fcs)^/ 

(fc + fc-2) + 

k2a n 



-9 

(-4) 



Act ng 

' fc+ 



fc+ 
p+ 



•"9 



(34) 



for connecting the same Feynman diagram in Fig. [TJb) to the second set tri-gluon correlation function Tag,f in 
Eq. The color factors in Eqs. ([55)1 and Eqs. ([M)) are the same as those in Eqs. and Eqs. (|32p . respectively. 

For our calculation of the evolution kernels in the next section in the light-cone gauge, n ■ A — the cut vertices 
are simplified as. 



fc+ 



p+ 



P^ 



V, 



LC 
G,F 



x{x + 0:2) {-gap) 5\x-—]x2 5\x2--^][ie 



P^ 



V: 



LC 



7+7^ 



A9,F — 2p-(- 



fc+ 



P^ 



■>,LC 



fc+ 



P^ 



a:(.T + 0:2) Uef 5 [x - — ] X25 [x2 - (-4) [-.9.^1 C^^^'^ 



P^ 



(35) 
(36) 
(37) 
(38) 



for the tri-gluon correlation functions, '7^,f, '^gf\ TAq,F, '^ag^f^ respectively. 
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B. Factorization and evolution equations 



In order to derive the evolution equations and evolution kernels from the definition of the twist-3 correlation 
functions, we need to compute the perturbative modification to these correlation functions caused by the quark-gluon 
interaction in QCD [3^. For example, we need to calculate the diagram in Fig. [2]for extracting the flavor non-singlet 
evolution kernel of the quark-gluon correlation function. 





FIG. 2: Feynman diagrams that contribute to the flavor non-singlet cliange of tlie twist-3 quark-gluon correlation function 
where fi, p and c are Lorentz and color indices of gluon field operators, respectively. The lower part of quark and gluon lines 
are contracted to the cut vertex that define the quark-gluon correlation function. 



We first evaluate the perturbative change to all spin-dependent correlation functions, Tq^p, Tq p , T/\q.P; and p, 
because they are defined in terms of hadronic matrix elements and represented by the Feynman diagrams with proper 
cut vertices. We follow the standard QCD coUinear factorization approach to factorize the perturbative change to 
these spin-dependent correlation functions into short-distance evolution kernels convoluted with corresponding gauge 
invariant long-distance matrix elements or the correlation functions [H, Hi], |4^. From Eqs. (fTTj) . (fT3| . ((201) . 
(P^ . we then derive the evolution equations for the two sets of twist-3 correlation functions that are responsible for 
generating the SSAs in QCD coUinear factorization approach. 

We start with the flavor non-singlet change to the quark-gluon correlation function, Tq^p, as represented by the 
diagrams in Fig. [5] with the cut vertex in Eq. For the twist-3 correlation functions relevant to the SSAs, we 

are interested in the difference of the diagrams in Fig. [2] with hadron spin st and that with —st- Only survival 
leading twist matrix element from the top of the diagram in Fig. ^a.) after taking the difference is the transversity 
distribution that does not contribute to the change of the quark-gluon correlation function of masslcss quark due to 
the symmetry of time-reversal or simply due to an odd number of gamma matrices in the spinor trace. Following 
the same derivation and steps presented in Refs. [4l|, |42|, we find that the twist-3 contribution from the diagram in 
Fi g. Ef a) can be combined with the leading contribution of the diagram in Fig. [H^b) due to color gauge invariance 
[sl llll| . Long-distance physics of the combined contribution from two diagrams in Fig.[2]could be expressed in terms of 
four twist-3 long-distance matrix elements or correlation functions, Tj^^'^}^, as defined in Ref. [Tol |. where superscripts 
V and A represent the vector and axial vector current, respectively, and subscripts, D and F, refer to the standard 
QCD covariant derivative and field strength, respectively. The correlation functions, Tp and Tp, correspond to our 
spin-dependent quark-gluon correlation functions, Tg^p and T^q^p, respectively, while the other two functions could 
be obtained by replacing the field strength operator F^'^ by the covariant derivative operator Da-. As explained in 
Ref. the two correlation functions with the covariant derivative operator do not contribute to the SSAs. 

We now provide a detailed derivation of the projection operator for extracting the flavor non-singlet evolution kernels 
or the short-distance contribution from both diag rams in Fig. [S] There are two sources of twist-3 or subleading power 
contribution from the diagram in Fig.[2{a) [8l.l4ll [4^. One is from the transverse momentum expansion of the parton 
momenta entering the bottom part of the diagram and the other is from the spinor trace decomposition when the 
bottom part of the diagram is contracted by 7 • n instead of the leading 7 • P [il], E^l . The part from the transverse 
momentum expansion contributes to the matrix elements with a covariant derivative, T^'^, which do not contribute 
to the SSAs as discussed in Ref. The other subleading contribution of the diagram in Fig.[2)Ja) due to the spinor 
decomposition could contribute. Although the matrix element of the subleading term from the spinor decomposition 
of the diagram in Fig. [D^a) has only two quark field operators, it can be expressed in terms of a matrix element of two 
quark fields and a gluon field by applying the equation of motion [4l[. That is, this part of subleading contribution 
from the diagram in Fig.l^^a) can be represented by the same diagram in Fig.[2Ub) except that the partonic Feynman 
diagrams in the bottom part of the diagram are given by the diagrams with the contact interaction (sl. [4^ . Therefore, 
we can derive the full flavor non-singlet evolution kernels from the diagram in Fig. [2{b) with the understanding that 
the bottom part of the diagram also includes those with the contact interaction or the special propagator [43 | . 
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We represent the perturbative change to Tq^p from the diagram in Fig. [UJb) as 



dTq^pix, X + Xi,[iF, St) = 



(27r)s 



Tr 



T''{p,P2,P,ST)Hp{p,P2,X,X2,fJ.F) 



(39) 



where the "Tr" represents the trace over the fermion fields' spinor indices, and T and H represent the top part and 
the bottom part of the Feynman diagram, respectively. In the momentum space, the T is given by the matrix element. 



TPip,p2, P, st) - (P, stI - P2)AP'%P2) ^g.jip) \P, st) , 



(40) 



where i,j are color indices of the quark fields and p,c are Lorentz and color indices of the gluon field, respectively. 
The H represents the bottom blob that includes all cut Feynman diagrams for the given external quark and gluon 
lines. The list of all cut diagrams at order of as will be given in the next section when we present the calculation of 
evolution kernels. The bottom quark and gluon lines of these diagrams are contracted by the cut vertex that defines 
the correlation function. The dependence of x and X2 in the argument of H in Eq. (j39p is from the cut vertex, and 
the scale fj^p represents the hardness or the off-shcUncss of the parton momenta, k and ^2. To pick up the leading 
power contribution from the perturbative modification to the quark-gluon correlation function, we first separate the 
spinor trace for the case of massless partons by [l^l 



1 



1 



Hp{j),P2,X,X2,HF) ~ Hp^a{p,P2,X,X2,^lF) yo^" j + H p,a{p, P2, X, X2, ^J,F) ( ^"f" [i"f 

where ". . ." represents terms with even number of 7-matrices and subleading, and 

1. 



Hp^a{p,P2,X,X2,IJ.F) = 2'^'' Hp{p,P2,X,X2,HF)lo 



Hp,a{P,P2,X,X2,fJ.F) = ^Tr Hp{p,p2,X,X2,fJ.Fha («7^ 



and 



In order to derive the contribution from the first term in Eq. (j41[) in details, we introduce 

d'^p d^p2 



T''" (p, P2 , P, ST)Hp,a {p,P2,X,X2,tJ-F) , 



(41) 



(42) 



(43) 



(44) 



with 



TP-''ip,P2,P,ST) = -TT 



TP{p,p2,P,STh' 



(45) 



We then apply the strong ordering in the off-shellncss of active partons, ^ and \p2\ ^ n'j?, and make the 
collinear approximation to expand the parton momenta entering into the H in Fig. [2l^b) around p = and p2 ~ 
as 



Hp^a{p,P2,X,X2,tJ.F) ~ Hp,a{^P,S,2P,X,X2,^J,F) 



dHp^a{£.P, ^2P,X,X2,flF) 



dpP 



{p-^Pf 



+ 



dHp^a{^P, ^2P, X, X2,Hf) 

dp2 



iP2 - £.2PY 



(46) 



By substituting Eq. (|46)) into Eq. (|44)) . we can rewrite the Iq as 



h - j d^d^2TP^''{t^ + a)HpA^,^2,x,X2,fiF)+ J d^d^2Tr'^{^,^ + ^2 



dHp,a{^,^2,X,X2, Hf) 

dpi^ 



J did^2Tr'^it^+i2) 



dHp^a{^,^2,X,X2,flF) 

dp2 



(47) 
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where the exphcit P dependence in Hp, a is suppressed. The correlation functions in Eq. (j47p are given by 

d-^pd% 



(27r)8 
d^pdH 



P2 



p+ 



7 



(2^) 



d'^pdH2 



(2^) 



Me 



P2 
P+ 



ip2 ~ ^2PY 



X (P, stI - P2)^AP'''{P2) i'qAp) \P, st) 



(48) 



Finally, we decouple the contraction of Lorentz indices in the RHS of Eq. (|47|) to express the quark-gluon correlation 
functions in terms of the Tq p-, defined in Eq. (p7|) . so that we can factorize the leading term of the RHS of Eq. ((47|) 
into a convolution of the 7^,f and corresponding evolution kernel. We find 



Co 



6 



(i f^Tpnn^^ pa 



(49) 



where the factorization scale dependence is suppressed and the fermionic color projection operator Cq is given by 

{Cqr^, = 2/{Nl-\W)J^. (50) 

with the quark and gluon color indices ij and c as labeled in Fig. [21(a), so that CqCq = \. In Eq. ((49|) . the quark-gluon 
correlation function p has the same definition as that of Tq p in Eq. (|12p . except the cut vertex in the square 
brackets is replaced by the cut vertex in the light-cone gauge in Eq. (|35p . The superscript "LC" indicates that this 
quark-gluon correlation function is calculated by using the light-cone gauge cut vertex instead of the full cut vertex. 
We find that the term proportional to Pf '"''^(^, ^ + ^2) hi Eq. (|T7)) does not give the leading power contribution. For 
the third term in Eq. (|47p . we have 



Ca ('ie"^''"") PP P" 



7;5°)(?,e + 6,5T) 



(51) 



where the long-distance quark-gluon correlation function, T^'f''^ ^ same definition as that of Tq p in Eq. (I12p . 

except the cut vertex in the square brackets is replaced by 



2P+ 



p- 



P+ 



P^ 



Pi \ P2a rip 



P2 



(52) 



which corresponds to the second term in the square brackets in Eq. ([29|) . The superscript "CO" indicates that this 
term provides the leading contribution in a covariant gauge calculation of the correlation functions d, HI] . From the 
factorized expression for the first and the third term, we find the leading contribution from the RHS of Eq. (|47p can 
be factorized as 



d^d^2i^if\£.,^ + i2.ST) 

+7;^^°He,C + 6,ST) 



-1 



C, ( — ) (ze^-''"«)P"ifp,„(f,6,2:,X2,MF) 



Cq (ie 



STf3 nn\ 



pP po 



dHp^a{.L£,2,X,X2,HF) 



dp2 



+ . . . , (53) 



where the ". . ." again represents the subleading term which includes the contribution from the Tf "''^ in Eq. g2 
From the definitions of T^^\^,^ + C2,st) and T^^\^,S. +^2,st), we have 



r.A^, ^ + 6, st) = r'^""' (e, ^ + 6, st) + rjT\^, S + 6, st) • 



(54) 
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Therefore. QCD color gauge invariance requires 



-1 
1^ 



STpnn\ pQ /^(LC) 



pp p 



9P2 



(55) 



when the LHS is evaluated in the light-cone gauge and the RHS is evaluated in a covariant gauge. Then, the two terms 
in Eq. (|53p can be combined into one term proportional to the quark-gluon correlation function, Tq^p{^,S, +S,2,st)- 
Since T^'^\S,,£, + ^2:St) vanishes in the light-cone gauge, the left-hand-side (LHS) of the equality in Eq. ([55]) 
represents the short-distance partonic part calculated in the light-cone gauge. On the other hand, the RHS of the 
equality in Eq. ([55]) represents the short-distance partonic part calculated in a covariant gauge [1] . This is because the 

matrix element T^'^\^,^ + ^2,st) dominates over T^^\^t ^ + ^2, st) in a covariant gauge calculation [l,!!!]. That 
is, the equality in Eq. (|55p provides an excellent consistency test for the perturbative modification of the quark-gluon 
correlation functions evaluated in different gauges. 

By using Eqs. (|54p and (j55p . we can combine the two factorized terms in Eq. (|53p into one factorized term as 



d£.2 %i,Fi£., £,+S.2, St) dKqq{i, i + £,2-,X,X + X2,fJ-F) 



(56) 



where the perturbative modification to the correlation function, dKqg{^,£^ + ^2,x,x + X2,fJ-F), can be calculated by 
using either side of the equality in Eq. ([55)1 depending on the gauge used for the calculation. For the light-cone gauge 
calculation. 



dKqq{^,£, + ^2,X,X + X2, I^f) = ^9 ( ^ 



From Eqs. (jH]) and Eq. ([57)) . we derive the projection operator in the light-cone gauge. 



V. 



(LC) 
<1,F 



:^l■P[^]{^e^-P-^)Cq, 



(57) 



(58) 



for extracting the perturbative modification dKqq from the partonic diagram in Fig. EKa) , which is equal to the lower 
blob of the diagram in Fig. Eljb) plus all diagrams with the contact interaction. From the RHS of Eq. (|55|) , we have 
the projection operator for the covariant gauge calculation 



(CO) 



-/■P pp (ie"^''"") Cq 



where the p2 is set to following the derivative 
^P>c 





d 



P.c 



(59) 



fk+ky 





FIG. 3: Partonic Feynman diagrams that contribute to the evolution kernels of the twist-3 correlation functions. 

In order to derive the leading contribution from the second term in Eq. (|4ip . we introduce 

d'^p d^p2 



Aq 



(27r)8 



TP'°'{p,P2,P, ST)Hp^aiP,P2, X, X2, Mf) ■ 



(60) 



with 



TP'''{p,P2,P,st) = ^Tr [fp{p,p2, P, st)i" (*7') 



(61) 
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Following the same derivation as that for /„, we find 



I 



pp p 



dp2 



d^2 TAq,Fi^,£,+ S.2,St) dKqAqid + ^2,X,X + 2:2, ^_f) 



where the perturbative modification to Tq.f from T^q^F is given by 

dKqAq{£,.i + S,2.X,X + X2,y^F) = Cq(^^^{i4)P''Hj,':^\^,^2,X,X2,flF) 

a ^Hj,^aH^,^2,X,X2,^J.F) 



(62) 
(63) 

(64) 



= C„ 



pp p 



d/2 



P2=i2P 



where the subscript "LC" ("CO") again indicates the light-cone (covariant) gauge calculation. From Eq. (|64p . we 
obtain the projection operator, 



-r.(LC) ^ n 5 
T^Aq.F ^1^1 -Pi 



6 



( S?p) Cq , 



(65) 



for extracting dRq^q from the same diagram in Fig. [Sj^a) in the light-cone gauge. Similarly, one can easily derive the 
projection operator for the covariant gauge calculation from Eq. (j64p . 

By adding contributions from Eqs. (|56p and ([63|) . we obtain the factorized perturbative modification to Tq^p, 



dTq^Fix, X + 2:2, ^_F, St) 



d^d^2 X},f{C,^ + ^2,St) dKqq +i2,X,X + X2,^J.F) 

+T'Aq,F{£,, C + ^2, St) dKqAqiC, £.+£,2,X,X + X2, ^If) 



(66) 



As shown in the next section, the leading power perturbative modification, dKij with i,j = q, Aq, g, Ag, can be 
expressed as 



i-Pf (lf^2 

dKij{^, ^ + ^2,X,X + X2, flF) = J -j^ C + ^2,X,X + X2,as) 



(67) 



where Kij {(,,^ + ^2,x,x + X2,as) is referred as the short-distance perturbative evolution kernel. Substituting Eq. (p7|) 
into Eq. (|66p and taking the derivative with respect to the factorization scale /ij? in both sides in Eq. (|66p . we derive 
the leading order flavor non-singlet evolution equation for the quark-gluon correlation function. 



^i%^^Tq,F{x, X + X2, /iF, St) 

oa4 



d^di2 



Tq,F{i, C + 6, Mi^, St) Kqq{^, ^ + £.2, X, X + X2, Us) 
T'\q.F{^, £. + ^2,fiF, St) KqAqii, £, + (,2, X , X + X2, Us) 



(68) 



which is consistent with that in Eq. ^ and has the generic homo gene ous differential-integral form of the typical 
evolution equation, such as the DGLAP evolution equation of PDFs |39l. |45| . 

Next, we derive the perturbative change of tri-gluon correlation function Tqp^ from the diagrams in Fig. ID Since 
gluon transversity distribution vanishes [49| . there is no leading twist or leading power contribution to the evolution 
of the tri-gluon correlation functions from Fig. IDJa). Similar to the case of the flavor non-singlet change of quark- 
gluon correlation function discussed above, the subleading power contribution from the diagram in Fig. [UJa) can be 
combined with the leading contribution of the diagram in Fig. Hl^b) [43 |. We can then derive the projection operator 
for calculating the gluonic evolution kernel by factorizing the diagram in Fig.lJJb). 

We express the diagram in Fig. [H^b) as 



dTQ^p{x, X + X2 



,f^F,ST)^ j '^-0^ [r'''"'^(p,P2,P,ST)i?^l^(p,P2,X,X2,MF3 

P,a,l3 



'AG ' 



(69) 
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FIG. 4: Feynman diagrams that contribute to the change of the twist-3 tri-gluon correlation functions where a, (3, ji, p and 
a, b, c axe Lorentz and color indices of gluon field operators, respectively. The lower part of gluon lines are contracted to the 
cut vertices that define the tri-gluon correlation functions. 



where the superscript i = f,d from the cut vertex and Iq (Iag) represents the part of the perturbative change that 
is symmetric (antisymmetric) in the exchange of the Lorentz indices a and /3. In Eq. (|69p . the partonic part Hp^a,p is 
given by the bottom part of the diagram in Fig. IH^b) plus diagrams with the contact interaction from the subleading 
contribution of the diagram in Fig. SJa). All partonic diagrams are contracted by the cut vertex in Eq. (|32p. The 
tri-gluon matrix element TP'°'-P in Eq. is defined as 



(70) 



with the gluon color indices, 6, c, a and is represented by the top part of the Feynman diagram in Fig. Hl^b). Following 
the same steps used to factorize the diagram in Fig.l^Jb), we can factorize the leading power contribution to the part 
that is symmetric in a and (3 as 



r(0 



^G.F (^,$ + $2,St) 

?(C + 6) 

^G,F Ui?+?2,Srj 



6 



Hp%^f3{S.,i2,X,X2,^J-F) 



(71) 



s 2 ^ ' dvl 



where the transverse polarization tensor d'^^ = —g"'^ 

1 



with j = f,d are given by 



cb 



1] 



if act, and (d''^)ach 



n^n^ and the gluonic color projection operators C, 



Xi) 
3 



(7V2-4)(7V2-1) 



^acb ■ 



(72) 



for color indices labeled in Fig.HJb), so that C^^^ C^^'^ = 1 for j = /, d. In Eq. 



^^f^) and fi^i(^°) with, = 

are tri-gluon correlation functions that have the same definition as that of T^p^ in Eq. pO|) . except that the cut 
vertex is now replaced by the corresponding one in the light-cone gauge and the one in a covariant gauge, respectively, 
and they satisfy 



^G:r(e ? + 6, 5t) - r^['^>'--> it e + 6, st) + r^^p 

Again, the color gauge invariance requires 



1 



-1 



dp^2 



(73) 



(74) 



P2=e2P 



when the LHS is evaluated in the light-cone gauge and the RHS in a covariant gauge. Therefore, the two leading 
terms in Eq. (|7ip can be combined together as 



l]^^ ^ j di d^2 T^%{t i + 6, St) dR<;p {t^ + t,X,X + X2,liF) 



with 



,STpnn\ 7:7-(i)(LC) 



(75) 



(76) 
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in the light-cone gauge. The dK^p can also be derived in a covariant gauge from the RHS of Eq. ([74]) . The equality 
in Eq. (|74p provides an independent check of perturbative calculation done in different gauges. From Eq. (j76p . we 
obtain the light-cone gauge projection operator, 



1 



-1 



^(/.d)(LC) _ 1 ,„/3 1 — I u e 



(77) 



for calculating the perturbative modification to the tri-gluon correlation function p from the diagrams in Fig.[3]^b). 
which includes all the partonic Feynman diagrams from the lower blob of the diagram in Fig. IH^b) plus corresponding 
twist-3 contribution from the diagram in Fig. SJa) expressed in terms of diagrams with the contact interaction [i^ . 
Similar projection operator can be derived from the RHS of Eq. ((74)) for the covariant gauge calculation. 



Similarly, we derive the perturbative modification to Tqp^ from the tri-gluon correlation function T^q\p, 



'AG 



iv) . 



(78) 



with 



(79) 



in the light-cone gauge. One can easily derive the expression for dK^^Pg in a covariant gauge as well. From Eq. (|79p . 
we obtain the light-cone gauge projection operator. 



7(/,'i)(LC) 
AG,F 



(80) 



for calculating the perturbative modification from the tri-gluon correlation function T^q\, from the same diagrams 
in Fig.[3);b). 

Using the generic expression of the leading power contribution to the perturbative modification factor dKij in 
Eq. (|67p ; we derive the evolution equation for the factorization scale dependence of the tri-gluon correlation function 

Tqp^ by factorizing the perturbative correction from the diagrams in Fig. 51 

IJ-F-^f^%{x,X + X2, fJ.F, St) = X! / '^^'^^2 f^^pi^,^ + ^2,tJ'F,ST)K^gl'\^,^ + ^2,X,X + X2,as) 



^ag,f(^: C + 6> MF, st) Kl'l'gi^, C + ^2,x,x + X2,a,) (81) 



where the superscript i, j = /, d, A'g^*^ and K^''^^ are evolution kernels that can be perturbatively calculated from the 
diagram in Fig. [3fb) with proper projection operators as discussed above. 





FIG. 5: Feynman diagrams that contribute to the change of the twist-3 quark-gluon correlation function where a, j3, ^, p and 
a, b, c are Lorentz and color indices of gluon field operators, respectively. The lower part of quark and gluon lines are contracted 
to the cut vertex that defines the quark-gluon correlation function. 



The evolution equation for the scale dependence of the quark-gluon correlation function in Eq. (j68p can also get 
contribution from the tri-gluon correlation functions via the diagrams in Fig. (5] Similarly, the evolution equation 
for the tri-gluon correlation functions in Eq. (jSip can get additional contribution from the quark-gluon correlation 
function via the diagrams in Fig. [51 
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FIG. 6; Feynman diagrams that contribute to the change of the twist-3 tri-gluon correlation functions from the interaction 
initiated from the quark-gluon correlation functions. The lower part of gluon lines are contracted to the cut vertices that define 
the tri-gluon correlation functions. 



Following the same procedure to factorize the diagrams in Fig. |3J we derive the additional contribution to the 
evolution of the quark-gluon correlation function from the tri-gluon correlation functions and have, 



2 d ~ 



d£.d£,2 Tq,F{£., i + 6: St) Kgg{^, ^ + ^2,X,X + X2,as) 

+TAq,F{L £. + 6, MF, St) KgAqiS,, ^ + (,2,X,X + X2,as) 



+ 



E 

i=f.d 



d^d^2 



TgW, ? + 6, MF, st) A'«(e, ^ + ^2.X,X + X2,a,) 



,x,x + X2,as)\ . (82) 



The evolution kernels from the tri-gluon correlation functions to the quark-gluon correlation function, /^^g''^'' and 

^qAg ' can be obtained by calculating the diagram in Fig. [3jc) with proper projection operators. If the kernels are 
evaluated in the light-cone gauge, the three gluon lines on the top of the diagram are contracted by the projection 
operator in Eqs. ([77)1 and (|80p . respectively. The diagram in Fig.[2)Jc) includes all Feynman diagrams from the bottom 
part of the diagram in Fig. \^h) plus diagrams from the sublcading contribution of the diagram in Fig. ^a,) , which 
can be effectively expressed in terms of the diagrams with the contact interaction and the same external lines as those 
in Fig.[5l^b). The combination of these diagrams forms a gauge invariant set [i^ . 

Similarly, following the same procedure to factorize the diagrams in Fig. [21 we derive the additional contribution 
to the evolution equation of the tri-gluon correlation functions from the quark-gluon correlation functions and have, 

fJ'F-g-^'^G^Fi.^^^ + ^2,tJ-F,ST) = E / '^^^'^2 fQ pi^,^ + ^2,^J.F,ST)K^gf\^,^ + ^2,X,X + X2,as) 
' j=f,d 

+ ^AG,F ' ? + 6 , MF , st) ifj Ag ^ + ^2, X, X + X2, a^) 



(83) 



+ E / ^^^^^2 \Tq^F{^.^ + ^.2,^JiF.ST)K'^^{^,^ + ^2.X,X+X2,as) 

q 

+^a<},f(C, C + 6 , Mf, St) K^^Ag i^,^ + ^2,x,x + X2,a 



where runs over all quark and antiquark flavors, the superscript, i,j = /, d. The evolution kernels from the quark- 
gluon correlation functions to the tri-gluon correlation functions, Kg^'"^^ and K^g^^g , can be obtained by calculating 
the diagram in Fig. [S^d) with proper projection operators. If the evolution kernels are evaluated in the light-cone 
gauge, the quark and gluon lines on the top of the diagram are contracted by the projection operator in Eqs. ([58[) 
and ([65[) , respectively. The diagram in Fig. [31(d) includes all partonic Feynman diagrams from the bottom part of the 
diagram in Fig. [6jb) plus the diagrams with the contact interaction representing the subleading contribution of the 
diagram in Fig. [H^a). 

Following the same derivation for the pcrturbativc corrections to the first set twist-3 correlation functions, we derive 
the evolution equations for the second set of twist-3 correlation functions, 



y?F-^^'^Aq,F{.X, X + X2, M_F, St) 

uHf 



d^d^2 



'?Ag,F(Ci C + 6: MF, St) KAqAqit ^ + ^2,X,X + X2,as) 
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+ ^AG,f (C, ^ + 6, Aif , St) K^AqAgii^ ^ + ^2,X,X + X2, , (84) 



and 



l4^g-^'^AG,Fi^^^ + ^2,fJ'F,ST) = ^ d£_ d£_2 T^q ,p(^, C + ^2, MF, St) i^AgAglC' + 6, 2: + X2, a^) 

+ ^G,Ut S + ^2, MF, st) K^Aggi^^ ^ + ^2,X,X + X2,a,) 
+ X! I d^di2 fq^F{£.,S.+£.2,^^F,ST)K'{!gq{^,S,+£,2,X,X + X2,as) 

q 

+T'Aq,FiL^ + ^2,t^F,ST)K''ll^g{^,^ + ^2,X,X + X2,as) . (85) 

All evolution kernels in Eqs. ([84]) and (|85|) can be derived by calculating diagrams in Fig. [3] with proper projection 
operators discussed in this section. 

Equations ([5^ . and ([55]) form a closed set of evolution equations for the scale dependence of the two sets 

of twist-3 quark-gluon and tri-gluon correlation functions defined in the last section. From these evolution equations, 
we can construct the evolution equations of twist-3 correlation functions that are responsible for the SSAs as, 



fIp^^Tq^F{x,X + X2,^_f) 



d 



^p-^^T^^p{x, X + X2,^iF) 
Ollp 

AiF^^^Ag,F(a;,a; + X2,tiF) 



d 



|J:i,—^T^%p{x,X + X2,^J-F) 



d 



d 



fip—-^Tq^F{x, a; + a;2, AiF, St) + fJ^p-^^Tq^p^x + X2,x, ^J,F,ST) 



d 



fJ-F-^^T'G%{x, X + X2, /iF, St) + UpT^T^'pix + X2,x, ^ip, St) 



d 



9/i 
d 



(86) 
(87) 



d 



fip^-^TAq,Fix,X + X2,^J-F,ST) - fJ.p^-^TAq,F{x + X2, X, flp, St) 



/^F^^^AG,f(2;, X + X2, ^iF, St) - A^F ^^^AG,F (^^ + X2,X, fXp, St) 



(89) 



As we show in the next section, the sum or the difference in the RHS of above equations determines the symmetry 
property of these correlation functions when the active momentum fractions a; and x + X2 are switched. 



IV. EVOLUTION KERNELS 



We present in this section our calculation of the order of evolution kernels for the evolution equations that are 
derived in the last section at X2 = 0. More precisely, we derive the order of as evolution equations for the diagonal 
twist-3 correlation functions defined in Eq. We will present the complete evolution kernels at the order of as in 

a future publication. 

The evolution kernels can be derived from the order of diagrams in Fig. [3] after setting 2:2 = or integrating over 
X2 weighted by 6{x2). We use the light-cone gauge cut vertices and projection operators derived in the last section to 
contract the quark and gluon lines at the bottom and the top of these diagrams, respectively. Since the cut vertices 
with the middle gluon in the LHS of the cut are the same as that with the gluon in the RHS of the cut, we only 
need to calculate the cut Feynman diagrams in Fig. [3] that have the middle gluon at the bottom part of the diagrams 
in one side of the cut. On the other hand, the sum of the all final-state cuts requires us to calculate all diagrams 
with the middle gluon on the top part of the diagrams in both sides of the cut. In addition, we need to calculate the 
same diagrams in Fig. [3] with the active momentum fractions x and x + X2 switched, as indicated by the equations in 
Eq. dMD to Eq. 1^. 

We start with a detailed calculation of the order of as evolution kernels for the evolution equations of 7q^p(x,x + 

X2, fJ-F, St) s^T^d %i,Fix + X2, X, fip, St), and then, we construct the evolution equation for 7^. F(a^, 2;, /if) from Eq. (|86p. 
Finally, from Eq. (|15p . we have the diagonal correlation function, Tg p(x,x,iip) = 2T:Tq p{x, x, fip). We define 

dlqq= dX2S{x2) dKqq{^,£, + £_2,X,X + X2, as) , (90) 



18 




FIG. 7: Feynman diagrams that contribute to the leading order flavor non-singlet evolution kernel of the twist-3 quark-gluon 
correlation function. 



where dKqq is given by the diagrams in Fig. [Sj^a) with the cut vertex in Eq. p5|) and the projection operator in 
Eq. (|58p . We list in Fig. [7] all cut Feynman diagrams at order of as with the gluon at the cut vertex in the LHS of 
the cut. Diagrams labeled from (a) to (m) have the top middle gluon in the LHS of the cut while the diagrams from 
(n) to (q) have the top middle gluon in the RHS of the cut. The quark propagator with a short bar for the diagrams 
labeled by (1), (m), (n), and (o) is the special propagator introduced in Ref. [43 to represent the contact interaction. 
These diagrams represent the contribution from the diagram in Fig. [DJa) that is necessary to make the full twist-3 
contribution gauge invariant. In the n ■ A = light-cone gauge, the Feynman rule for the special quark propagator of 
momentum k is 



42| 



17 • n 



2fc • n fc2 



(91) 



Having the cut vertex and the projection operator for the bottom and top quark and gluon lines, respectively, 
calculation of these Feynman diagrams in Fig. [7] is straightforward. In particular, after setting X2 = or integrating 
over X2 weighted by the S{x2), all diagrams labeled from (f) to (q) give the vanishing contribution to the diagonal 
evolution kernel, Kqq{^, ^ -t- ^2, x, x, as). Using the technique introduced in Ref. (soj . we find the following results for 
the rest of diagrams. 
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(93) 



(94) 



(95) 



In above equations, z ~ x/^ and the color factor for each diagram is explicitly shown in the square brackets with 
Cp = {N^ — l)/2Nc, Ca = and Nc ~ 3, the number of color. We notice that the RHS of the last equation from 
the diagrams (d) and (e) is infrared divergent, and the divergence is needed to cancel the infrared divergence from 
the term proportional to Cp in Eq. (|92[) when — > 1 . This cancellation of infrared divergence between the real and 
virtual diagrams is the same as that takes place in the evolution kernel of normal PDFs [39[ . The remaining infrared 
divergence as z — > 1 in Eq. (j92p is proportional to a different color factor, Ca/2, and is cancelled by the contribution 
from diagrams (b) and (c). 

From the same Feynman diagrams in Fig. [71 we can also calculate the contribution from T^q^p, 

(96) 



dlqAq = / dX2 6{x2)dKqAq{S,,^ + £,2,X,X + X2,as) 
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by using the projection operator in Eq. (j65p . In this case, only diagrams (b) and (e) give nonvanishing results, 

'Ca 



dl. 



(b) 
qAq 



2tt 



1 1 

2 x + £_2 
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Ois 
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27r 





(97) 
(98) 



By comparing above calculated results with Eq. (j67|) . we extract evolution kernels, Kqg{(^, S.~^S,2iX,x) and Kq/^q(^, ^ + 
^2tX,x). By calculating the same diagrams in Fig. [3] with momentum fractions ^ and x switched with ^ + ^2 and 
X + X2, respectively, we derive evolution kernels, Kqq{£^ + ^2, ^, x, x) and Kq/^qi^ + ^2,£,^ x, x). By integrating Eq. ([5^ 
over X2 weighted by 5{x2) or simply setting X2 = 0, we obtain the order of Ua evolution equation for Tq^F{x, x, fip) 
from flavor non-singlet interactions. 



dTq^F{x,x,y.F) 



2tt 



'^1 



Pqq{z)Tq,F{£„i,^iF) 

14 



Ca 
' 2 

Ca 
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[%,F{^,X,^J.F) - Tg,F(C,'f,MF)] + zTq^F(^,X,^iF) 



TAq,F{x,£,, f^p) 



(99) 



where 



Pqqiz) = Cf 



1 + ^' 3^. 



(100) 



is the LO quark-to-quark splitting function for the normal PDFs. The standard definition of "+" distribution is 



dz 



dz 



1 



+ /(1) In(l-x) 



(101) 



for a smooth function f{z). In deriving Eq. Eqs. (|lip and were used. It is clear from Eq. ([M)) that 

the flavor non-singlet evolution kernels for the diagonal twist-3 quark-gluon correlation function Tq^pix, x, ^f) = 
2TTTq,Fix, X, fip) are all infrared safe. The evolution equation for the diagonal correlation function Tq^p^x, x, fip) is 
not a closed one since it gets contribution not only from the same diagonal function Tq^F{£,,£,, fJ-p) but also from the 
off-diagonal part of the same function as well as gets the contribution from a different function T^q^pix, fxp). 

In the rest of this section, we derive the order of ttg evolution kernels involving gluons as well as those with the 
flavor change. In Fig. [SJ we list all cut Feynman diagrams at the order of as that could contribute to the evolution 



kernels, Kg^p and -fi^XgAg '^ith «,j = /, d, when proper cut vertices and projection operators are used. The gluon 
propagator with a short bar in the diagrams (1), (m), (n), and (o) is the gluonic special propagator defined in Ref. [i^ . 
which represents the contact interaction. The diagrams with the contact interaction are responsible for the twist-3 
contribution from the diagram in Fig. [UJa). We calculate all diagrams with the cut vertices and projection operators 
derived in this section and setting 0:2 = 0. We find that after taking 2:2 = or integrating over X2 weighted with 
6{x2), only diagrams (a), (b), (c), (d), (e), (f), and (g) give the nonvanishing contribution to the evolution kernel. 
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FIG. 8: Feynman diagrams that contribute to the leading order evolution kernel from the tri-gluon correlation functions to the 
tri-gluon correlation function. 



where Uf = 1, 2, ... the number of active quark flavors and the color factors arc shown in the square brackets. We find 
that the results from Eq. (|102|) to (|106p are the same for the evolution kernel K^^^^ and Kgg'^\ All evolution kernels 
of the crossing contribution Kgl"^^"^ = KfP = 0. 




FIG. 9: Feynman diagrams that contribute to the leading order evolution kernel from the tri-gluon correlation functions to the 
quark-gluon correlation function. 




■tp+P-, 




FIG. 10: Feynman diagrams that contribute to the leading order evolution kernel from the quark-gluon correlation functions 
to the tri-gluon correlation functions. 



By adding the flavor changing contribution to the evolution kernels from Figs.lOlandllOl and adding the contributions 
from the same diagrams but with their parton momentum fractions ^ and ^ + ^2 switched, we derive the order of as 
evolution equations for the diagonal twist-3 quark-gluon and tri-gluon correlation functions defined in Eqs. ([1]) and 
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In above evolution equations, z = a;/^, the LO quark-to-quark splitting function is given in Eq. (jlOOp . and the rest 
LO parton-to-parton splitting functions for the normal PDFs are given by 



= \ [(1-^)' + ^'] , 



P,q{z) = 2Ca 
PM = Cf 



1 - z 



(l-z)+ 

(1-Z)2 + 1- 



+ z(l - z) 



(Ill) 



Our explicit calculation also verifies that evolution equations for the diagonal parts of the second set of twist-3 
correlation functions vanish, 



dTAq(q),Fix,X,flF) 



= 0, 



dT^Q^p{x,x,HF) 



= 0. 



(112) 



din ' i91n 

which are consistent with the antisymmetric nature of the second set of twist-3 correlation functions. The off-diagonal 
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correlation functions in Eqs. (|107p to (|110p arc defined as 

Tq{q),Fi^,x,^iF) = 27rT,(^)^^(^,x,^i^) , rAg(g),_F(C,a;,Mi=^) = 27rTAg(g)^F(C,a;,/iF) , 

T^y{tx,^J.F) = 27r — ^ ^ , T^-'pj,(^,x,/ii.) = 27r ^ — , (113) 

where the more general correlation functions Tq^p(S,,x, ^p), T^GF^i^^XjjjLF), T^Aq,F{^-,x, (jlf), and 'T^^\,{S,,x, hf) are 
defined in Eqs. ([TT|) . (fT3|) . ((20|) . and ([22]) . respectively and the extra 27r is due to the fact that the / dj/^ in Eqs. ([1]) 
and docs not have l/27r. 

Equations from (|107p to (|110p plus Eq. (|112p form a complete set of evolution equations of the diagonal twist-3 
correlation functions relevant to the leading giuonic pole contribution to the SSAs. All evolution kernels at the order 
of Us are infrared safe and perturbatively calculated. However, unlike the evolution equations for the full twist-3 
correlation hmctions from Eq. ([55]) to Eq. ([55]) . these evolution equations do not form a closed equation set. 

The evolution equations of the diagonal twist-3 correlation functions have a lot in common with the DGLAP 
evolution equations of normal PDFs. Every channel of parton-to-parton evolution kernel is led by a term that is 
proportional to the DGLAP splitting function and the diagonal twist-3 correlation functions. All other terms are 
either proportional to the difference between the diagonal and the off-diagonal correlation functions or proportional 
to the off-diagonal correlation functions. Therefore, we expect that the scale dependence of the diagonal part of the 
twist-3 correlation functions is more close to the scale dependence of spin-averaged PDFs, not the spin-dependent 
helicity distributions. 

Unlike the normal PDFs, the quark-gluon and antiquark-gluon correlation functions could have a different evolution 
equation unless the tri-gluon correlation function Tq p = 0. The difference was caused by the difference in color 

contraction for Tq \ and T^f^p. As pointed out in Ref. [1^, the production of open charm mesons in semi- inclusive 
deep inelastic scattering or hadron-hadron collisions can provide the excellent information on the tri-gluon correlation 
functions. If Tq p ^ 0, the difference between the quark-gluon and antiquark-gluon correlation functions could be 
enhanced as the scale evolves. The difference should lead to interesting measurable consequences when we compare 
the SSAs generated by the quark-gluon correlation with that by the antiquark-gluon correlation. 

It was argued in Ref. ^2^ that one of the two tri-gluon correlation functions T^'^p can be related to the moment of 
a TMD gluon distribution, known as the giuonic Sivers function, in terms of their operator definitions. However, the 
other tri-gluon correlation function T^'^p does not have a direct operator connection to the TMD gluon distribution. 

Equation piOp indicates that within QCD coUinear factorization formalism, the Tq^p can be generated perturbatively 
from the quark-gluon and antiquark-gluon correlation functions as long as Tq^F{x,x, ^f) + Tq^pix, x, ^f) ^ or 
-^AG f(^'^'A*f) 7^ 0, even if T^^^p vanishes at one scale. 

To complete this subsection, we state that we also examined the infrared sensitivity of the order of as evolution 
kernels for correlation functions that give the leading fermionic pole contribution to the SSAs. The fermionic pole 
contribution is generated by the off-diagonal correlation functions, 7^,f(0, ^f), Th^FiOjX, ^f), ^ij,f(0, a;, 
and Tag. f{0,x, ^f) (or equivalently from 7^^i?(x, 0, /ii?), 7G,_F(a;, 0, /i^), T^q.Fix^O, ^f), SLXid T^g, Fix, 0, ^f) for the 
diagrams in which the gluon at the cut vertex is in the RHS of the cut). At the order of as, all evolution kernels are 
also infrared safe. For example, the flavor non-singlet evolution kernel for the quark-gluon correlation function Tq p 
can be calculated from the diagrams in Fig. [7] by setting a; = 0. We find that after setting a; = 0, only diagrams (a), 
(b), (c), (e), (f), and (g) in Fig. [7] give nonvanishing contribution to the evolution kernel. Again, the evolution kernel 
is infrared safe and all infrared divergences cancel between diagrams. We will present the full evolution equations for 
the off-diagonal twist-3 correlation functions in a future publication. 



V. SCALE DEPENDENCE 



In this section, we study the scale dependence of the diagonal twist-3 quark-gluon and tri-gluon correlation functions 
relevant to SSAs by solving the evolution equations derived in the last section. 

Since the evolution equations for the diagonal twist-3 correlation functions from Eq. (|107p to (jllOp do not form 
a closed set of differential equations, we need to make a model for off-diagonal correlation functions before we can 
study the scale dependence of the diagonal correlation functions. For the following numerical study, we introduce the 
following model to express the symmetric off-diagonal correlation functions in terms of diagonal correlation functions 
and a universal width, 

1 (xi-x2)^ 
Tq^F{xi,X2, ^J-F) = 7^ [Tq.Fixi,Xl,flF) +Tq^F{x2,X2,IJ,F)] e , 
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(114) 
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ofT-diagonal correlation. From Eq. pi3p . the ofF-diagonal correlation function Tqp\xi, X2, fJ-p) defined in Eq. (|113p 
is not symmetric. From Eq. (|114p . we have 



where both Tq^p and Tqp'^ are symmetric in exchange of xi and X2 and the cr is a width for the strength of the 



T, 



G,p {X1,X2, flp) = - 



We choose the width a such that 
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when l^i — X2I — > 1. In Fig. [Til wc plot the factor e 2^ as a function of a; for <t = 1/4 (solid line) and 1/8 (dashed 
line). 
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FIG. 11: The factor e 2^ as a function of x for cr = 1/4 (solid) and cr = 1/8 (dashed). 



To numerically solve the evolution equations in Eqs. p07p to piop . we choose the following input correlation 
functions at /ipo = 2 GeV. For the quark-gluon correlation function Tq^p{x,x, /j,fq), we choose the Fit. II of the 
quark-gluon correlation function Tq^p{x,x, ^p) from Ref. [1^. For the tri-gluon correlation functions, wc choose the 
model introduced in Ref. 12811. 



T^'p{x, X, upo) = A/ G{x, fipo) 



T^^p{x,x,fipo) = ArfG(x,/iFo) 



(117) 



with Xf = Xd ^ 0.07 GeV at /ii?o = 2 GcV and CTEQ6L unpolarized gluon distribution [50[. As an approximation, 
we also set T/\q^p(x,S,, iip) = and T/^q_p(x,S,, ^p) = for less parameters since they have vanishing diagonal 
contribution and the size of the off-diagonal part could be smaller than that of set one correlation functions. For a 
better convergence of the numerical solution, we use the linear combination of the two tri-gluon correlation functions, 



T, 
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Q p — p -1- p^ when we solve for the evolution equations. 
We first solve the flavor non-singlet evolution equation for the quark-gluon correlation function in Eq. ([M]) to test 
the relative role of the normal DGLAP evolution term that is proportional to Pqq{z) and the new piece that depends 
on the off-diagonal correlation function. In Fig. 1121 we plot the twist-3 up-quark-gluon correlation Tu^p{x, x, np) as 
a function of x at the factorization scale fip = 4 GeV (left) and fip = 10 GeV (right). The difference between the 
left figure and the one on the right indicates the evolution of the twist-3 correlation functions. The factorization scale 
dependence is a solution of the flavor non-singlet evolution equation in Eq. (|99p . Solid and dotted curves correspond 
to two different choices of the width for the off-diagonal input correlation function at cr = 1/4 and 1/8, respectively. 
The dashed curve is obtained by keeping only the DGLAP evolution kernel Pqq{z) when we solve the flavor non-singlet 
evolution equation in Eq. (|99p . Similarly, we plot the twist-3 down-quark-gluon correlation Td^p{x, x, fj.p) as a function 
of X at the factorization scale iip = A GeV (left) and ^p = 10 GeV (right) in Fig. [TSl Unlike the up-quark-gluon 
correlation function Tu,p, the down-quark-gluon correlation function T^^p is negative |8|. Figures [T^ and [T3] clearly 
show that the scale dependence of the diagonal twist-3 quark-gluon correlation function does follow the evolution of 
the unpolarized quark distribution. The difference between the solid and the dashed curves indicates that the effect 
of non-DGLAP type contribution from the off-diagonal correlation function could be very important at small x if the 
width of the off-diagonal correlation function is large. 
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FIG. 12: Twist-3 up-quark-gluon correlation Tu,f{x,x, fip) as a function of x at fip = 4 GeV (left) and fip ~ 10 GeV (right). 
The factorization scale dependence is a solution of the flavor non-singlet evolution equation in Eq. (|99p . Solid and dotted curves 
correspond to a = 1/4 and 1/8, while the dashed curve is obtained by keeping only the DGLAP evolution kernel Pqq{z) in 
Eq. 




FIG. 13: Twist-3 down-quark-gluon correlation Td^F{x,x, ^f) as a function of x at fip ~ 4 GeV (left) and ~ 10 GeV 
(right). Solid and dotted curves correspond to a = 1/4 and 1/8, while the dashed curve is obtained by keeping only the 
DGLAP evolution kernel Pqg{z) in Eq. (f99)l . 



In Figs. [T4landll5| wc plot the twist-3 up-quark-gluon and down-quark-gluon correlation functions. Tutrix, x, ij,f) 
and Td^F{x,x, hf), as a function of a; at /ij? = 4 GeV (left) and fiF = 10 GeV (right). Only difference between the 
solid and dotted curves in these figures and those in Figs. [12] and [13] is that we use the full set of evolution equations 
in Eq. p07|) through piOp to solve for the factorization scale dependence of these correlation functions. The dashed 
curves represent the quark-gluon correlation functions obtained from the parametrization of Fit II in Ref. [3^ by 
assuming all quark-gluon and tri-gluon correlation functions obey the DGLAP evolution. Wc find that non-DGLAP 
terms in the full evolution equations for the diagonal twist-3 correlation functions play a significant role in modifying 
the evolution of these correlation functions at small x, where the role of the off-diagonal correlation functions is 
enhanced due to a larger available phase space for the evolution kernels. The extra enhancement of the solid and 
dotted curves over the dashed curves in Figs. [TJ] and [15] is mainly from the term proportional to the sum of both 
tri-gluon correlation functions T^'^p and T^^p that we assumed to have the same sign. 

In Figs. [TBI and [TTl we plot the twist-3 tri-gluon correlation functions, Tqp{x,x,^f) and Tq\p{x,x, (Xf), as a 
function of a; at = 4 GeV (left) and ^f = W GeV (right), respectively. Solid and dotted curves are from solving 
the full evolution equations with the input correlation functions evaluated at cr = 1/4 and 1/8, respectively. Dashed 
curves are given by the normal CTEQ6L gluon distribution multiplied by the normalization constant Xf (or A^), 
which corresponds to making an assumption that all twist-3 correlation functions obey the DGLAP evolution, like 
the normal unpolarized PDFs. We notice that for the evolution of tri-gluon correlation functions, the difference in 
color factor for the DGLAP-typc terms in the full evolution equations tends to compensate the contribution from the 
terms proportional to the off-diagonal correlation functions, so that the evolution of the tri-gluon correlation functions 
follow more closely to the DGLAP evolution as shown in Figs. [TBI and [TTl 

We complete this section by stressing that the scale dependence presented in this section is sensitive to our as- 
sumption to neglect the role of the second set of twist-3 correlation functions and our model for the input tri-gluon 
correlation functions (equal and positive at the input scale) . Although the overall features found here should be valid. 
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FIG. 14: Twist-3 up-quark-gluon correlation Tu,f{x,x, ^.f) as a function of x at — 4 GeV (left) and fiF ~ 10 GeV (right). 
The factorization scale dependence is obtained by solving the full set of evolution equations in Eq. (|107p through (|110|l . Sohd 
and dotted curves correspond to a = 1/4 and 1/8 for the width of input off-diagonal correlation functions. The dashed 
curves represent the quark-gluon correlation functions obtained from the parametrization of Fit II in Ref. [33l ] by assuming all 
quark-gluon and tri-gluon correlation functions obey the DGLAP. 




FIG. 15: Twist-3 down-quark-gluon correlation Tu,f{x, x, (if) as a function of x at hf = 4 GeV (left) and /if ~ 10 GeV (right), 
all curves are defined in the same way as those in Fig. 1141 



the precise numerical values of these correlation functions should be extracted from a consistent global QCD analysis 
by comparing experimental data on SSAs and corresponding theoretical calculations, like what have been done to 
test the leading power QCD factorization formalism 0, 0- The new evolution equation derived in this paper is the 
necessary step to make such a consistent global QCD analysis possible for twist-3 correlation functions relevant to 
SSAs. 




10'^ 1 10'^ 1 

X X 



FIG. 16: Twist-3 tri-gluon correlation function Tq p{x,x, fiF) as a function of s at fiF ~ 4 GeV (left) and /if = 10 GeV (right). 
Dashed curves are from Tqp{x, x, /if) = A/ G{x, /j-f), and solid and dotted curves are from solving the full evolution equations 
with a — 1/4 and 1/8 for the input correlation functions, respectively. 
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FIG. 17: Twist-3 tri-gluon correlation lunction Tq^p{x,x, ^f) as a function of 2; at ~ 4 GeV (left) and fiF ~ 10 GeV (right). 
Dashed curves are from Tq''p{x,x, /if) ~ XdG{x, f^F), and solid and dotted curves are from solving the full evolution equations 
with a — 1/4 and 1/8 for the input correlation functions, respectively. 

VI. SUMMARY AND CONCLUSIONS 

We constructed two sets of twist-3 correlation functions that arc responsible for generating the SSAs in the QCD 
coUinear factorization approach. These correlation functions do not contribute to the long-distance correlation func- 
tions extracted from any measurement of parity conserving double-spin asymmetries, such as the DIS structure 
function (72- We introduced the Fcynman diagram representation for the twist-3 quark-gluon and tri-gluon correlation 
functions and derived the cutvertices to connect the hadronic matrix elements of these correlation functions to the 
forward scattering Feynman diagrams. In terms of the Feynman diagram representation, we derived for the first time 
a closed set of evolution equations in Eqs. ([M)) and ([55]) for these quark-gluon and tri-gluon correlation functions. We 
also provide the explicit prescription for calculating the corresponding evolution kernels. 

We calculated in the light-cone gauge the order of ag evolution kernels for the scale dependence of the diagonal 
quark-gluon and tri-gluon correlation functions that arc responsible for the leading gluonic pole contribution to the 
SSAs. We also provided in this paper the cut vertices and corresponding projection operators necessary for performing 
the calculation of evolution kernels in a covariant gauge. A covariant gauge calculation requires additional Feynman 
diagrams that have eikonal lines connecting quark and gluon lines at the cut vertices. These additional diagrams 
are from the expansion of the gauge links in the matrix element definition of the twist-3 correlation functions 0]. 
Our calculation explicitly shows that all evolution kernels are infrared safe as they should be. In addition, we also 
examined infrared behavior of the evolution kernels that are more relevant to the fermionic pole contribution to the 
SSAs and found that these kernels are also infrared safe. We leave the details on the scale dependence of more general 
off-diagonal twist-3 correlation functions to a future publication. 

We found that the evolution equations for the scale dependence of the diagonal twsit-3 correlation functions do 
not form a closed set of differential equations. The variation of the diagonal correlation functions gets contribution 
from the off-diagonal correlation functions. At the order of ag, all contributions to the evolution kernels from the 
diagonal correlation functions arc proportional to corresponding DGLAP evolution kernels of unpolarized PDFs. In 
order to test the role of the off-diagonal correlation functions to the scale dependence of the diagonal correlation 
functions, we introduced a model for the off-diagonal correlation functions in Eq. (|114|) in terms of the diagonal 
functions plus a correlation factor characterized by a parameter cr. With our model, we were able to solve the 
evolution equations and derive the scale dependence of the diagonal twist-3 correlation functions, and test the relative 
strength of different sources of contribution to the scale dependence. We found within a reasonable parameter space 
of our model that the scale dependence of the diagonal twist-3 quark-gluon and tri-gluon correlation functions follows 
closely the DGLAP evolution of unpolarized PDFs, and the contribution from the off-diagonal correlation functions 
becomes more significant at small x. 

Although our observations are dependent of the model for the off-diagonal correlation functions, we believe that 
the overall feature found here should be valid. The precise behavior of these twist-3 correlation functions should be 
extracted from a consistent global QCD analysis, just like what have been done for the leading power QCD dynamics 
and the scale dependence of the universal PDFs 0| . The evolution equations derived in this paper enable us to 
evaluate the NLO corrections to the SSAs systematically, which represents a necessary step moving toward the goal of 
global analysis of QCD dynamics beyond what have been explored by the very successful leading power QCD coUinear 
factorization formalism. With such a global analysis, we will be able to explore the rich field of quantum correlation 
of multiple fields inside a polarized hadron. 
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